Superconducting and Pseudogap effects on the interplane conductivity and Raman 
scattering cross section in the two dimensional Hubbard Model 
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Cluster dynamical mean field methods are used to calculate the superconductivity-induced changes 
in the interplane conductivity and Raman scattering cross section of the two dimensional Hubbard 
model. When superconductivity emerges from the pseudogap, the superconducting response is found 
to be diminished in amplitude, broadened and, in the case of the interplane conductivity, shifted 
to higher frequency. The results are in agreement with data on high temperature copper-oxide 
superconductors indicating that the Hubbard model contains the essential low energy physics of the 
pseudogap and its interplay with superconductivity in the cuprates. 
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I. OVERVIEW 



Three important characteristics of the layered copper 
oxide materials such as YBa2Cu306-i-2; are the existence 
of a large-gap insulating phase when a; = 0, of a "pseu- 
dogap" regime involving a suppression of the density of 
states, for a range of a; > 0, and of a d2;2_y 2 -symmetry 
superconducting state occurring for a range of x which 
partly overlaps the pseudogap regime. Understanding 
the relation (if any) between these phenomena is one of 
the central issues in the field. Two experimental probes 
which have been important in the discussion of this issue 
are the interplane (c-axis) conductivity and the Big Ra- 
man scattering cross section. These spectroscopies are 
interesting because they appear to be controlled by the 
"antinodal" electrons which are most affected by super- 
conductivity and the pseudogap and exhibit dramatic 
variations with temperature and carrier concentration. 
Both of these probes reveal striking temperature and 
x-dependcnt effects associated with the pseudogap and 
superconductivity, which one would like to understand 
theoretically. The problem has attracted considerable 
attention, but work to date has been primarily based on 
approximate analytical approaches, often involving par- 
ticular phenomenologically chosen ansatze for the rele- 
vant physical processes. 

The development of cluster dynamical mean field 
theorjiir— has changed the theoretical situation, provid- 
ing an unbiased (in the sense of not pre-selecting a par- 
ticular interaction channel or set of diagrams) numeri- 
cal approach to determining the properties of some of 
the basic models of condensed matter physics. Recent 
algorithmic developments^"— enable implementations of 
this procedure which are now at the point of provid- 
ing a semi-quantitative solution for the normal-state^ 
and superconductingi^ properties the two-dimensional 
square-lattice Hubbard model, which is widely believed 
to capture the essential physics of the high transition- 
temperature superconductors. Crucially, the new meth- 
ods enable access to clusters large enough to provide con- 



fidence that the true behavior of the Hubbard model is 
revealed, and to low enough temperatures that the su- 
perconducting state can be constructed. 

Dynamical mean field methods have determined that 
the Hubbard model exhibits a "pseudogap"i^iii"— where 
some regions of the Brillouin zone are gapped and oth- 
ers are noli^"— as well as a dj.2_j,2-symmetry supercon- 
ducting stated ' ^"i*^"" —. The superconducting state can 
now be constructed and aspects of its interplay with the 
pseudogap can be studiediS. Quasiparticle properties'^ 
and energetics^ have been determined. In this paper we 
use the methods to study the superconductivity-induced 
changes in the interplane conductivity and Big Raman 
scattering cross section of the two-dimensional Hubbard 
model. 

The rest of the paper is organized as follows. Sec- 
tion |TT] presents the response functions to be computed, 
gives more specifics of the physical phenomena of inter- 
est and outlines the theoretical methods used. Section 
mil presents the main new results of the paper: a com- 
putation of the interplane conductivity and Raman scat- 
tering cross sections. Section ITVl analyses the results and 
their relation to experiment. Section |V] summarizes the 
findings and outlines future directions for research. Ap- 
pendices provide calculational details. 



II. INTRODUCTION 

A. Model 

The essential structural motif of the high transition 
temperature copper oxide superconductors is the Cu02 
plane, a square planar array of Cu ions, with an oxygen 
ion at the midpoint of each Cu-Cu bond. It is by now 
accepted that the interplane coupling is weak enough that 
it may for most purposes be neglected, so that the basic 
physics problem which must be understood concerns the 
motion of electrons in a two dimensional lattice with a 
square symmetry. The interplane conductivity may then 



be studied by second order perturbation theory in the 
interplane couphng. 

In the "parent compounds" of high- Tc superconduc- 
tors, the density of electrons is one per Cu02 unit, but 
the materials are insulating with a large (^ 1.5 — 2 eV) 
band gap^. The insulating behavior is widely supposed 
to be a consequence of strong electronic correlations, re- 
lated to the "Mott insulating" phenomenon^. Remov- 
ing electrons (adding holes) produces metallic and su- 
perconducting behavior. The superconducting state is 
of d^2_y2 symmetry^, with the superconducting gap be- 
ing maximal at the center of the zone face ( (tt, 0) ) and 
vanishing along the zone diagonal ((0,0) — ?► (tt, tt)) di- 
rection. Adding holes also produces a region of "pseudo- 
gapped" behavior—"—, in which even at temperatures 
above the superconducting transition temperature the 
electronic density of states is suppressed in the zone face 
region but not in the zone diagonal region. 

A widely, but not universally, accepted hypothesis^ 
is that the basic theoretical model which describes the 
physics of interest is the two dimensional one-orbital 
square lattice Hubbard model, an idealized description 
which arises as a low (w < 1.5eV) energy effective model 
of the underlying material Hamiltonian. This model rep- 
resents the physical situation in terms of electrons mov- 
ing among sites of a two dimensional square lattice, and 
subject to an interaction U which disfavors double oc- 
cupancy. In a mixed position/momentum representation 
we have 



H = 



ka 



£kCl„Cka + 



uj: 



ni^nn- 



(1) 



J 



Here q^. creates an electron of spin a ~^, l on site i of 
a two dimensional square lattice of unit lattice constant, 
^k<7 '^^ ^^^ Fourier transform to momentum space, Ek is the 
energy dispersion and rii-f = cLci-|- is the operator den- 
sity of up-spin electrons on site i. In the computations 
presented below we take e^ == —2t (cos k^ + cos ky) 
with t = 0.35eV. While the dispersion is particle-hole 
symmetric we consider non-zero dopings for which the 
particle-hole symmetry is broken. 

Because we wish to treat superconducting phenomena 
it will be convenient to write subsequent equations in 
terms of the Nambu spinors 



and anomalous components of the electron Green func- 
tion, Eq. [31 We obtain these using the dynamical clus- 
ter approximation (DCA)ii^ version of cluster dynami- 
cal mean field theory. In this method the Brillouin zone 
is divided into some number TV of equal volume patches 
labeled by a central momentum K and the electron self 
energy T,{k,uj) is represented as a piecewise constant 
function taking different values in each patch, so defin- 
ing 4)K{k) = 1 if fc is in the patch labeled by K and 
(t>K{k) — otherwise. 



S(fc,w) 



K 



.K(fc)SK(w). 



(4) 



The YiK are obtained from the solution of an N-site quan- 
tum impurity model specified by the original Hubbard 
interaction and a self-consistency equation which relates 
lattice {G{k,uj)) and impurity {Gk{'^)) quantities: 



Gk{^)^N 



_fk_ 
{2ny 



cj)K{k)G{k,u). 



(5) 



The solution of the impurity model in the normal 
and superconducting states is obtained using continuous- 
time quantum Monte-Carlo methods^ with submatrix 
updates^ii as described in detail in Appendix \^ These 
methods provide solutions on the imaginary time or Mat- 
subara frequency axis; real frequency information is ob- 
tained from maximum entropy analytical continuation as 
described in Appendix [B] 

The computational burden of DCA calculations rises 
rapidly with cluster size A'' and interaction strength U . 
With the computational resources available to us, study 
of temperatures below the superconducting transition 
temperature at generic dopings is feasible for C/ < 7 at 
N = 8 but to obtain data of the precision needed for 
analytical continuation of response functions we employ 
U ~ 6. Fig. [1] plots the transition temperature against 
carrier concentration for this U. The onset of the nor- 
mal state pseudogap is also indicated as a dashed line. 
As will be discussed in the conclusions this U probably 
is slightly lower than the U which is relevant to the real 
materials. In consequence the superconducting region is 
pushed to lower carrier concentrations than observed in 
the real highTc materials. 



K = 



C-kl 



(2) 



and the corresponding matrix Nambu Green function de- 
fined for imaginary time r > as 

G(fc,r) = -(vI't(r)vI/,(0)\. (3) 



B. Formalism: dynamical mean field method 

Evaluation of the interplane conductivity and Raman 
scattering amplitude require knowledge of the normal 



C. Raman Big and Raman B2g cross section 

Raman scattering is a process in which an incident 
photon at some frequency cdin polarized along some di- 
rection a is scattered to an outgoing photon at some other 
frequency ujout = ^in — ^ and some other polarization di- 
rection b. The details of the Raman scattering process 
are complicated. Both phonons and electrons may be im- 
portant. The different choices of incident and outgoing 
polarization lead to many different symmetry channels. 
Often, to enhance the signal, the incident or outgoing 
photon frequencies are tuned to be near resonance with 
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FIG. 1. Solid lines and filled circles (black online): supercon- 
ducting critical temperature of the Hubbard model with near- 
est neighbor hopping calculated for U = 6t using the 8-site 
dynamical cluster approximation (see also Ref . VUX) . Dashed 
line (red online) denotes crossover to normal state pseudogap. 
Dotted lines (blue online) indicate the temperatures studied 
in this paper. Note that the temperature axis is cut off slightly 
below the lowest temperature used in this study. 



challenging. With the methods used here it has been re- 
cently attempted for a wide frequency range analysis of 
the normal state Raman cross section^. The method is 
too expensive to be run in the superconducting state on 
today's computers but the results of Ref. |40| indicate that 
while vertex corrections have a crucial effect at low dop- 
ings and high frequencies {'^ O.beV where they express 
e.g. the two-magnon contribution) they are much less 
important at the low frequencies of interest here. We 
therefore neglect the vertex corrections, so that on the 
Matsubara axis and in Nambu notation 

X(fi„) = T ^ Tr [TZlr^Gik, w„ + f]„)T3G(fc, w„)] . 

(9) 
The presence of the T3 factors in the expression for the 
Raman reflects the fact that the Raman probe does cou- 
ple to a quasiparticle at the gap edge. We evaluate the 
response function on the imaginary axis and then ana- 
lytically continue the result (see appendix [Bj) . 



D. Interplane conductivity 



some other excitation of the solid. Providing a detailed 
treatment of all of these effects is beyond the scope of 
this paper. Here we will focus on electronic scattering 
in the "-Big" channel at moderate energy transfers {il 
small compared to bandwidths or interband transition 
energies). In this case the basic Raman process is the 
creation of a particle-hole pair of negligible net momen- 
tum. For Big symmetry the amplitude for this process 
vanishes for electrons along the zone diagonal and is max- 
imal for electrons at the zone face. We will also present 
a few results for i?2g symmetry, where the amplitude is 
maximal for zone-diagonal electrons and vanishes at the 
antinodes. We follow the convention in the literature and 
assume the simplest forms compatible with the desired 
symmetries, which are 



i?i 



Raman 



-E^ 



k'^k,(y'^ko 



kcr 



with (in Big geometry) 



and (in i?2g geometry) 



TZ {coskx ~ cos ky) 



= TZsinkx sinky. 



(6) 



(7) 



(8) 



The amplitude TZ is not important for our purposes. 
Eas. l7l8l mav be derived in the non- resonant limit by use 
of a minimal coupling ansatz and are widely used in the 
literature. 

Standard linear response methods may then be used 
to obtain the susceptibility x(^) which characterizes the 
Raman response. The result involves both the electron 
Green function and a vertex function describing the in- 
teraction of the particlc-holc pairs created by the Ra- 
man process. Computation of the vertex function is very 



The computation of the interplane coupling begins 
from an ansatz for the coupling between planes. Elec- 
tronic structure calculations indicate that the coupling 
is small)^"— so only electron hopping between adjacent 
planes need be considered. This is parametrized by a 
hopping amplitude t±{kx,ky) which connects an elec- 
tronic state at in-plane momentum (k^, ky) in one plane 
with an electronic state of the same in-plane momen- 
tum in an adjacent plane. (Some authors have con- 
sidered impurity-mediated interplane coupling^ which 
does not conserve momentum, but we will restrict at- 
tention to the ideal undisordered situation here). Gauge 
invariance considerations then imply that the coupling 
of electrons to an electric field in the interplane direc- 
tion may be determined by multiplying the interplane 
hopping by the usual Peierls phase factor involving the 
vector potential A = Acz with Cz denoting the direc- 
tion perpendicular to the planes. Band theory and tight 
binding^— considerations indicate that t±_ has a strong 
momentum dependence for simple "one-layer" systems 
such as La2_a;Cu04 or TlBa2Cu06 so the final result is 
(setting e = h = c =interplane distance= 1) 



t±{k^ 



, ky] A) ~ —t±_ (cos kx — cos ky) e . 



(10) 



Thus the interplane hopping amplitude vanishes for mo- 
menta along the diagonals of the two dimensional Bril- 
louin zone and is maximal at the zone faces implying, as 
many authors have noted^^— , that the c-axis conductiv- 
ity is in effect a spectroscopy of the behavior of the zone- 
face electrons. In bilayer systems such as YBa2Cu306+i: 
the inter-bilayer coupling has the form given by Eq. [10] 
but the intra-bilayer coupling may also have contribu- 
tions from the zone diagonal electrons^. These compli- 
cations will not be addressed in this paper. 



In Nambu notation the Hamiltonian giving the inter- 
plane coupling is then 

(11) 

Here T3 is the Pauli matrix acting in Nambu space and for 
simplicity we have assumed ij_ to be the same between 
all planes. 

We now use standard linear response theory to write 
an expression for the conductivity valid to leading non- 
trivial order in both the applied vector potential and the 
interplane coupling. The conductivity is as usual the sum 
of paramagnetic and diamagnetic terms. In Matsubara 
space we have 



The c-axis supcrfluid stiffness pc^s is given bys^ 



with 



CTc(ao = f^f(ao + anao 



T 



(12) 



^f (^n) = Trll^r [ix(fc)V3G(fc,w™)r3G(fc,( 



i Lji 



rnk 



(13) 



T 



<^c{^n) = — ^Tr[ix(fc)'G(fc,c^™ + f]„)G(fc,w,„)] 



7nk 



(14) 

The absence of ra factors in the expression for a^ encodes 
the fact that a quasiparticle with energy equal to the su- 
perconducting gap is an equal admixture of electron and 
hole states, which because it is chargeless does not couple 
to an applied electric field. In materials with a bilayer 
structure extra contributions appear in the conductivity 
related to intrabilayer plasmon excitations^; these will 
not be considered here. 

In one important respect the piecewise constant self 
energy, an essential feature of the DCA dynamical mean 
field approximation, is problematic. As can be seen from 
Eq.[TU]there is an interplay between the magnitude of the 
superconducting gap ^ \kx — ky\ and the interplane cou- 
pling. Near the nodes the superconducting gap is larger 
than the interplane coupling, meaning that the contri- 
bution of nodal quasiparticles to the low frequency re- 
sponse is strongly suppressed. On the other hand, in the 
DCA approximation the entire momentum sector con- 
taining the nodal point is gapless, so that the Drudc re- 
sponse of nodal quasiparticles contribute strongly to the 
low frequency response. In our calculations we there- 
fore suppress the nodal region completely, by integrating 
only over momenta corresponding to the antinodal sec- 
tor. (This issue also arises in the Raman case, but is 
not important for the analysis because the response is 
suppressed at low frequencies). One may view our ap- 
proximation as using an interlayer coupling which has 
the same momentum discretization as the DCA approx- 
imation. 



,s^TY,t^{kfTi 



n.k 



T3G{k,UJn)T3G{k,UJn) 



G{k,UJn)G{k,UJn] 



(15) 



and the interplane conductivity may be rewritten as 



with 






T 



(16) 



a'.iQn) = ^f (^n) - o E ^^ (^)'Tr [G(fc, L.„).G(fc, c^„)] 



n,k 



(17) 

To obtain real-frequency conductivities we construct 
xd^n) = ^nO'ci^n), which is then analytically contin- 
ued using the methods employed for the Raman case as 
discussed in Appendix [Bland from this construct the con- 
tinued CTj, and thus ac- 

We remark that the interplane spectral weight 



Sc 



duj 



crc(w) 



(18) 



is given by i7„(T^(i7„)^ (note that the integral includes 
the superfluid delta function if present). 



III. RESULTS 
A. Raman Big and Raman B2g cross section 

Figure [2] presents one of the principal results of this 
paper: the doping dependence of the Big-symmetry Ra- 
man response of the two dimensional Hubbard model 
for several carrier concentrations at interaction strength 
U = 6i. This sequence of carrier concentrations corre- 
sponds to the cuts across the phase diagram shown as 
dotted lines in Fig. [1] The cuts include both Fermi liq- 
uid and pseudogapped regimes, and temperatures both 
above T = i/30 « imK (using t = 0.35eF) and well be- 
low the highest superconducting transition temperatures 
(T ^ i/60 w 7QK « r™'^^/2). 

The upper left panel of Figure [2] presents results ob- 
tained for a high doping just outside the regime where 
superconductivity is found at the temperatures we have 
studied. The high temperature Raman scattering am- 
plitude (solid curve, black online) has the features ex- 
pected of a strongly correlated Fermi liquid: a low fre- 
quency peak (visible at w « 250 cm^^) is characteristic 
of coherently propagating quasiparticles, while the rela- 
tively featureless higher frequency scattering intensity is 
attributable to the incoherent part of the electron spec- 
tral function. The peak frequency is related to the quasi- 
particle scattering rate. Indeed for quasiparticles with 
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FIG. 2. B\g Raman scattering cross section plotted against frequency (expressed in wave numbers using t = 0.35eV^) for the 
two dimensional Hubbard model with [/ = 6t at carrier concentrations 0.148 (top left), 0.102 (top middle), 0.075 (top right), 
0.053 (bottom left), 0.032 (bottom middle) and 0.015 (bottom right). Solid curve (black online) T = t/30; dot-dashed curve 
(red online) T = t/60, normal state; dashed curve (blue online) T — t/60 superconducting state. Label shading (color online) 
corresponds to color of curves for this doping in other figures. 



a momentum non-conserving scattering rate Tqp/2 the 
Raman cross section is 



xtm-xo- 



i + ( ji^ 



(19) 



As the temperature is decreased the low frequency peak 
moves to lower energy, indicating a decrease in scattering 
rate. For further discussion of the normal state issues see 
Refs. [23 andl40l and note that especially at lower dopings 
the frequency dependence at frequencies > 1500cm~^ is 
affected by vertex corrections not considered hero^. 

The upper middle and upper right panels of Figure [5] 
show dopings within the superconducting regime. For 
a; = 0.1 the normal-state Fermi liquid coherence effects 
are almost invisible. At the higher temperature T = t/30 
the scattering rate is large enough that the quasiparticle 
peak has merged with the continuum. As the temper- 
ature is lowered in the normal state a modest increase 
in the low frequency slope (decrease in scattering rate) 
is evident and a hint of a quasiparticle peak is seen as 
a very weak minimum ai Q, ^ 250c?7i^^, but there is no 
clear signature of coherent Fermi liquid behavior down to 
our lowest accessible temperature, x = 0.075 is a carrier 
concentration at the boundary of the pseudogap regime 
and near the maximum in the Tc vs doping plot. The 
normal state traces show a larger scattering rate (lower 
slope at r2 — > 0) and no trace of quasiparticle coherence 
is visible. 

At these carrier concentrations the onset of supercon- 
ductivity has dramatic effects. The low frequency in- 
tensity is suppressed, and a large peak becomes evident. 



This peak arises from quasiparticle excitations across the 
superconducting gap. It is visible in this response func- 
tion because the coherence factors are such that a quasi- 
particle at the gap edge couples to the Raman probe. In 
an s-wave superconductor with a uniform single-particle 
gap A the Raman intensity would diverge ^ l/\/fi — 2A. 
In a d-wave superconductor the momentum dependence 
of the gap eliminates the divergence, leaving just a peak, 
but in our calculation the piecewise constant nature of 
the self energy means that a divergence similar to the 
s-wave case will occur. However, the relatively strong 
inelastic scattering leads to substantial thermal broaden- 
ing of the divergence at the temperatures accessible to us. 
The peak structure still provides a good estimate of the 
gap 2A. At higher frequencies the superconducting state 
Raman amplitude appears to be smaller than the normal 
state amplitude, but these differences are at the edge of 
what can reasonably be resolved with presently available 
analytical continuation techniques so it is not clear how 
much significance can be attributed to this difference. 

The association between the peak in the Raman scat- 
tering intensity and the superconducting gap may be ver- 
ified from an examination of the electron spectral func- 
tion (imaginary part of the electron Green function di- 
vided by tt) shown for two carrier concentrations in Fig.|3l 
There is not a one-to-one correspondence, since the Green 
function is averaged over the (0, tt) momentum sector 
while the Raman intensity is in essence the average of 
a product of two G over the sector, and also involves 
coherence factors, but one may expect energy scales to 
be similar. The lighter weight lines (maroon online) show 
the normal (dashed line) and superconducting (solid line) 
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FIG. 3. Imaginary part of (0, n) sector electron Green func- 
tion (divided by tt) in normal (dashed lines) and supercon- 
ducting (solid lines) states for doping x = 0.102 (heavy lines, 
maroon online) and x — 0.032 (lighter lines, blue online). 
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FIG. 4. Big Raman scattering cross section in the normal 
(dashed) and superconducting (solid) state, for dopings x = 
0.075 (light lines, red online), x = 0.053 (intermediate lines, 
green online) and x — 0.032 (heavy lines, blue online). All 
curves are normalized to the value at 2000 cm~^ . 



spectral functions for the doping x ~ 0.102 whose Raman 
intensity is displayed in the upper middle panel of Fig. [5] 
No normal state pseudogap is visible. In the supercon- 
ducting state a clear gap is evident. The peak-to-peak 
separation ~ 600 cm~^ is seen to coincide with the peak 
position in the Raman intensity. 

The left and middle panels of the lower row of Fig. [2] 
display results for lower carrier concentrations, now well 
within the pseudogap regime. The pseudogap can be 
identified from the decrease in low frequency normal- 
state Raman intensity as temperature is lowered and 
a normal state pseudogap scale can be approximately 
identified from the frequency at which the lower tem- 
perature normal state trace rejoins the high temperature 
trace {^ 1100 cto~^ for x = 0.053 and ^ 1400 cm'^ for 
X = 0.032). These estimates are in reasonable accord 
with the pseudogap scales inferred from the normal state 
Green function. For example, the heavier dashed lines 
(blue online) in Fig. [3] show the normal state Green func- 
tion for doping x = 0.032. The pseudogap is clearly vis- 
ible in the normal state Green function and defining the 
pseudogap energy scale 2Apg as the peak to peak separa- 
tion in the spectral function gives 2Apg « 1400 cm~^ in 
agreement with the estimate from the Raman spectrum. 

For these dopings superconductivity again produces a 
peak in the Raman cross section, but the peak is broader, 
and the relative increase in amplitude less, than at higher 
dopings. To demonstrate this more clearly, in Fig. 2] 
we present the normal and superconducting state Big 
Raman scattering intensity, normalized to the value at 
2000 cm'^ for x = 0.102 (intermediate shading, ma- 
roon online) 0.053 (light shading, red online) and 0.032 
(dark shading, blue online). This rescaled plot makes it 
clear that as parameters are tuned through the pseudo- 
gap phase to the low doping superconducting boundary 
the change in Raman intensity due to superconductivity 
weakens both in relative and in absolute terms. 
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FIG. 5. B-ig Raman scattering cross section in the normal 
state at high T (solid line, black online) and at low T (dashed 
line, red online) and in the superconducting state at low T 
(dot-dashed, blue online), for carrier concentration x — 0.102 
(top panel) and 0.053 (bottom panel). 



We have also computed the Raman scattering cross 
section in the B2g scattering channel. The B2g matrix 
element (Eq. [5]) is maximal in the zone diagonal sec- 
tor, and contributions from the antinodal sector which 
carries the information about pseudogap and supercon- 
ductivity are small. Representative results are displayed 
in Fig. [S] The evolution of the spectra as temperature is 
lowered in the normal state is consistent with that re- 
ported previously2^. At all dopings a clear quasiparticle 
peak is observed which becomes better defined as tem- 
perature is lowered. The difference with the normal state 
Big spectra reflects the momentum-space differentiation 
characteristio^ii^ of the approach to the Mott transition 
in this model. The zone diagonal states which dominate 
the i?2g response remain more or less Fermi liquid-like 



while more exotic physics affects the states near the zone 
face which determine the Big response. The effects of su- 
perconductivity are difficult to discern in the calculated 
spectrum because in the DCA approximation used here 
the piecewise constant nature of the self energy means 
that the anomalous self energy vanishes in the entire 
zone-diagonal momentum sector where the B2g matrix 
element is peaked. A DCA calculation of the effects of 
superconductivity on the i?2g spectra would require a 
much finer momentum space resolution, corresponding 
to iV = 32 or larger, to capture the behavior of the su- 
perconducting gap near the nodes. These cluster sizes 
are not cannot yet be studied in the relevant range of T 
and U. 



B. Interplane conductivity 

Fig. ini shows the interplane conductivity calculated for 
the same parameters as the Raman scattering shown in 
Fig. [21 The normal state physics has previously been 
discussed in the context of 8-site DMFT calculations^^ 
and similar results for A^ = 2 and A^ = 4 have also been 
presentedji^— Our normal state results are consistent 
with these previous works. The overdoped case (upper 
left panel) exhibits a well defined low frequency "Drude" 
peak which rapidly sharpens as the temperature is de- 
creased. As the doping is decreased first the Drude peak 
is broadened and decreased in amplitude, then the tem- 
perature dependence ceases and a hint of pseudogap be- 
comes visible (upper right panel). With further decrease 
of doping the pseudogap becomes obvious; the pseudogap 
magnitude increases as the doping decreases. Note the 
large changes in y-axis scale between the two first panels, 
the panels for intermediate doping, and the panels for low 
doping, reflecting the drastic reduction of low-frequency 
interplane conductivity as doping is decreased although 
for frequencies greater than > 3000cm~^ the calculated 
interplane conductivity has only a weak doping depen- 
dence. 

Inspection of Fig. [6] shows that in the overdoped {x > 
0.1) regime, the low frequency spectral weight (integral 
of the conductivity over the region of the Drude peak) 
increases as T is decreased. On the other hand, in the 
underdoped regime {x < 0.07) the formation of the pseu- 
dogap corresponds to a decrease in the low frequency 
spectral weight. To determine if the spectral weight is 
drawn from other frequencies wc present in Fig. [7] the ki- 
netic energy, i.e. the integral of the optical conductivity 
over all frequencies. We see that for the normal state 
calculations, as temperature is decreased the kinetic en- 
ergy indeed increases on the overdoped side while on the 
underdoped side it decreases. We attribute the effect to 
temperature dependent modifications of the scattering 
rate for the antinodal carriers. 

We now turn to the effects of superconductivity. At 
X = 0.102 where the superconductivity emerges from a 
more or less Fermi liquid state, the opening of the super- 



conducting gap causes a suppression of the conductivity 
at low frequency (top middle panel). An increase in con- 
ductivity at a higher frequency ~ 600 cm^^ is also evi- 
dent. As the doping is decreased to the edge of the pseu- 
dogap regime {x = 0.075) and beyond, we see that the 
peak in the superconducting conductivity moves rapidly 
to higher energy, becoming comparable to the pseudogap 
energy scale. 

The existence of a superconductivity-induced peak 
is interesting because the BCS coherence factors^ are 
such that the conductivity vanishes at the supercon- 
ducting gap edge. Thus in standard BCS/dirty limit 
calculations of the c-axis optical response of a layered 
superconductor— this increase does not occur: at 57 f^ 
the superconducting state conductivity lies below the cor- 
responding normal-state curve. The changes therefore 
should be interpreted as arising from superconductivity- 
induced changes in the electron scattering rate, presum- 
ably associated with the pseudogap. Comparison of 
Fig.mto Fig.[2]shows that at the higher carrier concentra- 
tion the frequency scale in the peak in the superconduct- 
ing state c-axis conductivity matches that in the Raman 
cross section, while at the lower dopings the effects in the 
c-axis conductivity clearly occur at the pseudogap scale, 
not the superconducting gap scale. 

An important characterization of the superfluid prop- 
erties is the c-axis superfluid stiffness, shown for a range 
of dopings in Fig. [5] A clear maximum is visible. For 
dopings higher than the maximum, the well deflned 
Drude peak observed in the normal state means that con- 
ventional clean-limit (or intermediate-scattering) BCS 
physics applies. The decrease in penetration depth with 
increasing doping in this overdoped regime occurs be- 
cause our lowest accessible temperature is not far enough 
below the actual transition temperature, so that thermal 
excitations have reduced the superfluid stiffness. We be- 
lieve that if the calculation could be performed at very 
low temperatures, the superfluid stiffness would mono- 
tonically increase with doping until the high-doping end 
of the superconducting range. 

For the points x = 0.117,0.102 and 0.075 the acces- 
sible temperatures are sufficiently far below Tc that the 
results reflect the low temperature limit. It is important 
to note that for the value U = 6t studied here, x = 0.06 is 
approximately the point of maximum transition temper- 
ature and that all of these points lie outside the pseudo- 
gap regime. The decrease in penetration depth is there- 
fore not directly related to the pseudogap, but is due in- 
stead to the rapid decrease with decreasing doping of the 
low frequency normal state conductivity. As discussed 
in Ref. [2^ this decrease is due to the rapid increase of 
the low frequency scattering rate for antinodal electrons 
which is a precursor to the pseudogap. 

Comparison of Fig. [5] and [7] shows that for the dop- 
ings 0.12 > X > 0.06 where the low temperature limit 
is reached, the "FerrcU-Glover-Tinkham" sum rule is vi- 
olated, at about the 10% level. This means that the 
total optical integral in the superconducting state (in- 
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FIG. 6. Real part of interplane conductivity plotted against frequency (expressed in wave numbers using t = 0.35ey) for the 
two dimensional Hubbard model with [/ = 6f at carrier concentrations 0.148 (top left), 0.102 (top right), 0.075 (middle left), 
0.053 (middle right), 0.032 (bottom left) and 0.015 (bottom right). Solid curve (black online) T = i/30; dot-dashed curve (red 
online) T = t/60, normal state; dashed curve (blue online) T = t/&Q superconducting state (n > 0.02). Label colors (online) 
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eluding the delta function contribution) is different than 
that in the normal state. The sign of the violation de- 
pends on doping. On the overdoped side we see that 
the c-axis kinetic energy of the superconducting state is 
less than that in the normal state at the same temper- 
ature. On the other hand, in the underdoped state it 
is slightly greater. We attribute these effects mainly to 
the superconductivity-induced changes in the antinodal 
scattering rate (on the overdoped side) and changes in 
the pseudogap (on the underdoped side). 




FIG. 8. Superfluid stiffness pa determined in the supercon- 
ducting state at T = f/60 from Eq. ll5l as a function of doping. 



IV. DISCUSSION 

We begin the discussion by comparing our results to 
exp erimental datai^"— (for reviews see also Refs. l6a and 
67). As has been previously observed (see Refs. l25| and 
40D the evolution of the normal-state Big Raman ampli- 
tude with doping and temperature is in good agreement 
with data. Overdoped materials exhibit Raman spectra 
in reasonable accord with the spectra shown in the up- 
per left panel of Fig. [21 with a weak but visible coherent 
quasiparticle part which steepens as T is reduced. For a 
doping dependent temperature less than about 1507^' a 
quasiparticle peak appears, centered at r^ 200 cm~^. As 
the doping is decreased, first the peak vanishes as in the 
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FIG. 9. Difference of Raman intensity computed in the su- 
perconducting state at temperature T = t/60 and in the nor- 
mal state at temperature T — i/30 for carrier concentrations 
X = 0.075, X = 0.053, and x = 0.032. 



upper right panel of Fig. [2] and then the pseudogap leads 
to a suppression of the intensity over a wide range. See 
e.g. Fig. 1 of Ref. |6l. 

The superconductivity-induced changes are similarly 
in reasonable qualitative agreement with the data. For 
overdoped materials the onset of superconductivity leads 
to a large amplitude peak in the Raman intensity, at a 
frequency which it is reasonable to interpret as twice the 
gap magnitude. As doping is decreased, this peak moves 
to higher energy, loses intensity and broadens, and at 
very low doping the peak ceases to be visible. To make 
this more evident we plot in Fig. |9] the difference between 
the superconducting Raman spectra at our lowest tem- 
perature T = t/60, and the normal state Raman spectra 
at T = i/30. These trends are in reasonable qualitative 
agreement with the results presented in Fig. 4 of Ref. |6^. 

The theoretical curves shown in Fig. [2] reveal an addi- 
tional remarkable result: at the lower dopings the size of 
the superconducting gap (position of the maximum in the 
superconducting state Raman cross section) is less than 
the size of the pseudogap (energy scale at which the two 
normal state traces merge). This phenomenon may also 
be observed by comparing the peak-to-peak distances in 
the superconducting and normal state x = 0.02 spectral 
functions shown in Fig. [31 That the onset of supercon- 
ductivity leads to a decrease in the gap scale suggests 
that superconductivity and the pseudogap are compet- 
ing phenomena. The correspondence of the calculations 
to data suggest that this phenomenon also occurs in the 
actual materials. This idea is more extensively discussed 
elsewhereiSi^. 

One set of quantitative differences is that the numeri- 
cal values for gaps, for the dopings delineating different 
regimes, and the temperature scales are somewhat dif- 
ferent in the theory than in the experimental data. We 
believe this is a consequence of the relatively small value 
oi U = 6t which is numerically accessible to us. The re- 



sults of Ref. |lO| indicate that increasing U will increase 
the doping scales and decrease the gap values. However, 
Refs. |69J and |70| also report a coherence peak which re- 
mains reasonably sharp even in underdoped materials^, 
in contrast to the broadening observed here, and also 
indicate that, apart from the Big coherence peak, the 
temperature-dependent changes are of smaller magnitude 
than found in the calculations reported here. These are 
important issues for future investigation. 

We now turn to the interplane conductivity^ Here com- 
parison with experiments^— (see also Ref. 1821 ) is compli- 
cated because in "single layer" compounds the coupling 
between the planes is so weak that it is difficult to ob- 
tain reliable data while in the YBCO family of materials 
where the interplane conductivity is an order of magni- 
tude larger, bilayer plasmon effects complicate the in- 
terpretation of the data. Nevertheless several important 
points of comparison are possible. First, the qualitative 
feature that the conductivity is significantly doping de- 
pendent only at frequencies below ~ 2000cm~^ is consis- 
tent with data (see, e.g. Fig. 2 of Ref. |48() . Second, in 
optimal and overdoped materials the superconductivity- 
induced changes correspond to a strong decrease in the 
absorption at frequencies below ~ lOOOcm"^ with only 
a weak increase at higher frequencies > lOOOcm^^. As 
the doping is decreased the amplitude of the changes due 
to superconductivity (both the low frequency decrease 
and the high frequency increase) rapidly become smaller. 
However, an important difference is that the weak higher 
frequency peak does not shift as much with doping in the 
data as in the calculation. 

The rapid drop in superfiuid stiffness as doping is de- 
creased is consistent with observation (see, e.g. Fig. 2a of 
Ref. 1831) . Our calculations reveal that the initial stages of 
the drop are not due to the pseudogap per se, but instead 
reflect the dramatic increase in low frequency zone-face 
scattering rate which is a precursor to the pseudgap. 



V. CONCLUSIONS 

The high T^. copper oxide superconductors exhibit both 
dx2-y2 superconductivity, and a normal state "pseudo- 
gap" . In some regions of phase space these phenomena 
coexist and manifestations in different spectroscopies of 
the interplay between them has been of long-standing in- 
terest in condensed matter physics. The results of this 
paper, taken in conjunction with a wide range of previ- 
ous work)^"— strongly suggest that the two dimensional 
Hubbard model does contain the essence of the pseudo- 
gap and d-wavc superconductivity phenomena observed 
in the cuprates. The essential ingredient in the calcu- 
lation is the electron Green function, which is affected 
both by the pseudogap and by superconductivity. The 
interplay between these two phenomena, in combination 
with the BCS coherence factors, leads to the somewhat 
different behavior in Raman and c-axis conductivities. In 
the c-axis conductivity the key doping dependent changes 
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reflect a precursor of the pseudogap, namely the rapid in- 
crease in the scattering rate for electrons near the zone 
face. The correspondence between the calculated and 
measured Raman spectra lends support to the proposaliS 
that when superconductivity emerges from the pseudo- 
gap regime, the gap is decreased. 

The calculations presented here are not in precise cor- 
respondence to data. Because of limits on computational 
resources they are performed for an interaction strength 
U = 6t and with no second neighbor hopping t'. The 
onset of the pseudogap and the maximum in Tc occur 
at doping x w 0.07 rather than the w 0.14 observed 
experimentally^^. Available evidenc o^°i^^ indicates that 
the phase boundaries move to larger x as U is increased, 
but increasing U or t dramatically increases the severity 
of the fermion sign problem which is the crucial limit- 
ing factor in the calculation. Similarly, we studied the 
N = 8 cluster dynamical mean field approximation be- 
cause for larger cluster sizes studies at the necessary low 
temperatures and strong interactions are computation- 
ally too expensive. For N = 8 (and for N = 16, not 
studied here) the piecewise constant self energy used in 
the DCA dynamical mean field method produces a gap 
with three values: -l-A in the region of one antinode, 
—A at the other antinode and along the zone diagonal 
The much larger (iV = 25, 36) clusters needed to resolve 
the details of the momentum dependence of the gap are 
at present simply not accessible at the low temperatures 
needed for these studies. However, results presented so 
far in the literature^ make it clear that the N = 8, U = 6 
case studied here represents many essential features of 
the model. 

In summary, the physical picture emerging from these 
and other calculations is that the various charge-related 
experimental spectroscopies of the pseudogap and super- 
conductivity may be understood in terms of the pseudo- 
gap and superconductivity effects on the electron propa- 
gator, and that these effects may reasonably be studied 
in terms of the two dimensional Hubbard model. The 
important open theoretical issues are to extend the for- 
malism to the computation of magnetic quantities (which 
requires a treatment of vertex functions) and to under- 
stand the physical origin, in the model, of the supercon- 
ductivity and the pseudogap. More generally, relating 
these and related results to the more refined picture of 
the cuprates now becoming available, in particular to the 
growing evidence for the important of charge and spin 
stripe ordering^^— is an important open question. 
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Appendix A: Numerical Methodology 

We use a momentum-space (DCA) formulation of 
dynamical mean field theory, written following Maier 
et al.— in Nambu space, thereby allowing for supercon- 
ducting order—. We restrict our solution to the para- 
magnetic phase, suppressing long-ranged antiferromag- 
netism but allowing for antiferromagnetic fluctuations 
(albeit coarse-grained to the momentum resolution of the 
cluster). We focus on clusters of size eight; a size that is 
a compromise between accuracy (DCA becomes exact as 
N ^f oo but is approximate for any finite N) and numeri- 
cal expense (away from half filling, the cost of simulations 
as a function of U , N , and T^^ increases exponentially 
due to the fermionic sign problem). The N ~ 8 cluster 
approximation was found in previous work^ in the normal 
state to be large enough to distinguish generic iV — > oo 
behavior from that specific to particular clusters. Work 
in the superconducting state^^ showed similar behavior 
but large differences to simulations on smaller clusters 
(both CDMFT and DCA). 

The formalism requires an impurity solver formu- 
lated in Nambu space. Due to the small energy differ- 
ence between the superconducting and the normal state 
solution^ an unbiased, numerically exact solution of the 
impurity model is required. Continuous-time quantum 
Monte Carlo methods^i^ are the only impurity solver 
methods able to access to couplings strong enough to 
produce a pseudogap at temperatures low enough to 
construct the superconducting state numerically exactly. 
A variant of the Continuous-time Auxiliary field (CT- 
AUX) impurity solver with 'submatrix update' numeri- 
cal techniques^ makes such simulations feasible in prac- 
tice. The solver requires a decoupling of the (repulsive) 
Hubbard interaction in Nambu space, which is imple- 
mented analogously to the one in normal state presented 
in Ref. [y. The extension to Nambu space means that all 
matrices are twice as large as in a normal state compu- 
tation at the same temperature. 

We have found that the most stable procedure to ob- 
tain a converged solution is to begin at a relatively high 
temperature (e.g. T = t/10) and introduce a pairing 
field rji{k) = rjicj)!; via the replacement G(fc,ia;„;77i) = 
[za;„To -I- (^-efc)T3 -|-77i(fc)ri - S]" , with e.g. (/)fc = 
cos kx — cos ky for d-wave superconductivity, and rji typ- 
ically O.lt. Retaining the pairing field we obtain con- 
verged solutions G(-ftr, ia;„; 771) first at the initial temper- 
ature, then, using the solution at the initial temperature 
as a seed, at the desired range of lower temperatures. We 
remark that the sign problem for large 771 is less severe 
than at 77 = 0, so these computations are not inordinately 
expensive. 

Then, at each temperature, using the converged 
G{K,iuJn',r]i) as a seed, we set 771 = in the self- 
consistency condition and continue iterating until con- 
vergence is reached. At selected points we check the so- 
lution by taking the putatively converged self energy, di- 
viding the anomalous part by a large number (typically 
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20), and veriiying that under further iterations the solu- 
tion converges back to the one previously found (see also 
supplementary material of Ref. llOf ). 

The formalism produces normal and anomalous Mat- 
subara Green's functions G^/7v(iw„) and self-energies 
E^/^(ia;„), with constant relative errors of E: ^^ ^ 
const as a function of Matsubara frequency. 

We then use S^ and E^r to evaluate the expressions for 
the conductivity and the stiffness (Eq. [T7] and Eq. [T5)) 
and the 'bubble' term of the Raman response (Eq. |9]) on 
the Matsubara axis. We neglect the Raman vertex cor- 
rections that have been introduced in Ref. |40, based on 
the fact that they are (a) found to be small for small 
(real) frequencies, where we see the large superconduct- 
ing response, and (b) at the moment too expensive to 
compute numerically to the accuracy required for ana- 
lytic continuation. 



Appendix B: Analytical Continuation 

In the condensed matter physics context the analytical 
continuation problem is the inversion of a relation of the 
general form 



Miifln) 



dx S{x) 

TT in„ — 



(Bl) 



where M is a quantity measured on the imaginary fre- 
quency (or time) axis with measurement uncertainties 
which are relatively small, assumed to be Gaussian and 
characterized by a covariance matrix. 5* is the corre- 
sponding real axis spectral function. Because the ker- 
nel l/(i51„ — x) has many small eigenvalues, its inver- 
sion is an ill-posed problem. With the statistical meth- 
ods commonly used in condensed matter physics it is not 
easy to quantify the uncertainties in 5* arising from the 
combination of the approximate inversion and the mea- 
surement uncertainties in M. We view the process as 
one of data fitting, which generates a spectral function 
that is consistent with the Matsubara data within error 
bars. We find that if the measurement uncertainties are 
sufficiently small (our relative error is typically smaller 
than lO"**) different implementations of the continua- 
tion process produce a reasonably robust and consistent 
representation of data. In general the lower frequency 
structures are most reliable and small differences over 
large frequency ranges are less robust to variations due 
to choices in the continuation procedure. 

In the most widely studied case, Af is a measurement 
of the electron Green function on the fermion Matsubara 



points ujn — (2n + 1)ttT and S is the electron spectral 
function, which is non-negative and normalized so that 



dxS{x) — 1. 



(B2) 



We find that our covariance matrix for M(r2„) is almost 
diagonal when measured in frequency space,— so that 
correlations between different bins can be neglected. 

To continue the Raman scattering amplitude we 
observe that the Raman susceptibility also obeys a 
Kramers-Kronig relation 



X(*^n) = 



dx Imx(a;) 



(B3) 



but here the spectral function is odd: lrax{x) = 
— Imx(— x). To deal with this we reformulate the prob- 
lem as follows: 



X{i^n) 



dx Imx(a;) x 

■K X ifln — X 

dx Imx(a;) 



-1 



iflr. 



VQUf-t JU 



(B4) 

(B5) 



dx C(x) 
Xii^n = 0) + if^„ / — .^^ ' (B6) 



TT iQn — X 



with 



Cix) 



Imx(a;) 



We then invert the equation 

X(J^„) - xi'i^^i = 0) f dx C{x) 



iflr, 



n in„ 



(B7) 



(B8) 



using standard methods to find C, which is normalized 
according to 



dxC'ix) = xii^n = 0) 



(B9) 



and then construct the imaginary part of the Raman re- 
sponse as XRaman = wC(w). 

The c-axis conductivity is continued in a similar man- 
ner, defining Xad^n) = ^nC^(f^n) and then following 
the steps that led to Eq. IB6I 

We use the open source maximum entropy analytic 
continuation program available as part of ALPS^ for all 
of our continuations. 
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